1 Introducción

Este estudio presenta un análisis básico de RNA-Seq con énfasis en las posibilidades de visualización de datos y resultados^.

En su versión actual los gráficos son correctos pero poco elaborados. Para mejorarlos se utilizaran gráficos con ggplot2 utilizando como modelo algunas de las propuestas del documento: https://github.com/UofABioinformaticsHub/DataVisualisaton_BIS2016.git

1.1 Infraestructura informática para el análisis

Este estudio se llevará a cabo usando R/Bioconductor por lo que es preciso tener instaladas un conjunto de librerías. Esto puede hacerse siguiendo el procedimiento descrito a continuación.

El código que se presenta debería ejecutarse tan sólo una vez!

if (!require(BiocManager)) {
    install.packages("BiocManager", dep = TRUE)
}

installifnot <- function(pckgName, BioC = TRUE) {
    if (BioC) {
        if (!require(pckgName, character.only = TRUE)) {
            BiocManager::install(pckgName)
        }
    } else {
        if (!require(pckgName, character.only = TRUE)) {
            install.packages(pckgName, dep = TRUE)
        }
    }
}
installifnot("limma")
installifnot("edgeR")
installifnot("org.Hs.eg.db")
installifnot("clusterProfiler")
installifnot("dplyr", BioC = FALSE)
installifnot("gplots", BioC = FALSE)
installifnot("ggvenn", BioC = FALSE)
installifnot("pheatmap")
installifnot("prettydoc", BioC = FALSE)
installifnot("ggnewscale", BioC = FALSE)
# Pels grafics que empren ggplot2
installifnot("locfit", BioC = FALSE)
installifnot("magrittr", BioC = FALSE)
installifnot("statmod", BioC = FALSE)
# Pel HeatMap Interactiu
installifnot("heatmaply")
# Pel VolcanoPlot interactiu
installifnot("here")
installifnot("janitor")  # Cleaning column names  
installifnot("scales")  # Transform axis scales   
installifnot("ggrepel")

Puede resultar útil, aunque no es para nada imprescindible, disponer de, al menos dos, directorios específicos: - datos (o un nombre similar) donde guardar y de donde cargar los datos. - resultso (o un nombre similar) donde escribir los resultados.

if (!dir.exists("datos")) dir.create("datos")
if (!dir.exists("results")) dir.create("results")
if (!dir.exists("figures")) dir.create("figures")

2 Los datos para el análisis

Este análisis se basa en un estudio publicado recientemente (Arunachalam 2020) que investigaba la respuesta inmune a la infección con SARS-COV-2 desde un a perspectiva de biología de sistemas utilizando tecnología de secuenciación, RNA-Seq.

Los datos se han depositado en el repositorio “Gene Expression Omnibus (GEO)” con el identificador GSE152418.

2.1 Lectura de la matriz de contajes

Los datos de contajes se han descargado desde l repositorio “GEO” al archivo Rawcounts.csv que se utilizará como punto de partida para el análisis.

counts <- read.csv("datos/RawCounts.csv", row.names = 1)
selectedCounts <- as.matrix(counts)

En la misma web se dispone de información los grupos y otras covariables o variables auxiliares. Con éstos se crea un objeto (habitualmente denominado “targets”) que debe estar sincronizado con el anterior, es decir, cuyas filas deben corresponderse con las columnas de la matriz de datos.

sampleNames <- colnames(selectedCounts)
grupos <- c(rep("COVID", 17), rep("SANO", 17))
colores = c(rep("red", 17), rep("blue", 17))
selectedTargets <- data.frame(samples = sampleNames, group = grupos, cols = colores)
rownames(selectedTargets) <- selectedTargets[, 1]

3 Preprocesado de los datos

3.1 Estandarización de los contajes

Los datos de secuenciación pueden estar “desbalanceados” en el sentido que distintas muestras pueden contener un número distinto de secuencias, lo que puede inducir erróneamente a pensar que un gen se expresa más en una muestra que en otra, cuando esto se deba a esta diferencia global.

Esto puede evitarse expresando los contajes como “CPMs” es decir “counts per million”, lo que no modificará los resultados de comparaciones posteriores, pero hará que las muestras sean comparables en número , lo que es útil y necesario para los análisis posteriores.

library(edgeR)
selectedCounts[1:5, 1:6]
                COV145 COV147 COV149 COV150 COV151 COV152
ENSG00000223972      0      0      0      0      0      0
ENSG00000227232      1      0      3     16      1      2
ENSG00000278267      3      1      2      0      3      0
ENSG00000243485      2      0      0      1      0      0
ENSG00000284332      0      0      0      0      0      0
apply(selectedCounts, 2, sum)
  COV145   COV147   COV149   COV150   COV151   COV152   COV153   COV154 
15048983 13946337 13520020 18863644 13659353 13279390 12967503 14300187 
  COV155   COV156   COV157   COV175   COV176   COV177   COV178   COV179 
17192277 15626012 17041245 13415720 13679816 16092084 17379155 13672179 
  COV180   HEA061   HEA062   HEA063   HEA064   HEA065   HEA066   HEA067 
13722892 17003281 17359694 15570080 16172199 19420885 18264376 16398893 
  HEA068   HEA069   HEA070   HEA071   HEA181   HEA182   HEA183   HEA184 
14895860 14939005 20458479 17820834 14973734 14280768 18297901 17177433 
  HEA185   HEA186 
16467873 16860317 
counts.CPM <- cpm(selectedCounts)
counts.CPM[1:5, 1:6]
                    COV145     COV147    COV149     COV150     COV151    COV152
ENSG00000223972 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991 0.1506093
ENSG00000278267 0.19934902 0.07170342 0.1479288 0.00000000 0.21962973 0.0000000
ENSG00000243485 0.13289935 0.00000000 0.0000000 0.05301203 0.00000000 0.0000000
ENSG00000284332 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
apply(counts.CPM, 2, sum)
COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155 COV156 COV157 
 1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06 
COV175 COV176 COV177 COV178 COV179 COV180 HEA061 HEA062 HEA063 HEA064 HEA065 
 1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06 
HEA066 HEA067 HEA068 HEA069 HEA070 HEA071 HEA181 HEA182 HEA183 HEA184 HEA185 
 1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06  1e+06 
HEA186 
 1e+06

Una vez los datos estan como CPMs, se procede a filtrarlos,

3.2 Filtraje de genes poco expresados

Los genes con recuentos muy bajos en todas las librerías (es decir en todas las muestras) proporcionan poca evidencia de expresión diferencial, por lo que es habitual eliminar aquellos genes que, o bien varían muy poco entre grupos, o bien presentan poca o nula expresión en la mayoría de las muestras.

En este caso, siguiendo un criterio habitual, se opta por conservar únicamente aquellos genes que presentan algún valor en, al menos, tres muestras de cada grupo.

thresh <- counts.CPM > 0
keep <- (rowSums(thresh[, 1:10]) >= 3) & (rowSums(thresh[, 11:20]) >= 3)
counts.keep <- counts.CPM[keep, ]
dim(counts.CPM)
[1] 60683    34
dim(counts.keep)
[1] 33039    34

Aunque no sea más que un ejemplo basta con ver los dos primeros genes para comprobar como el primero no cumple la condición, en el grupo “SANO”, mientras que los siguientes sí que la cumplen. Por lo tanto, al filtrar desaparece el primer gen de la matriz filtrada, pero no los dos siguientes.

round(head(counts.CPM), 3)
                COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155
ENSG00000223972  0.000  0.000  0.000  0.000  0.000  0.000  0.000   0.00  0.116
ENSG00000227232  0.066  0.000  0.222  0.848  0.073  0.151  1.388   0.42  1.163
                COV156 COV157 COV175 COV176 COV177 COV178 COV179 COV180 HEA061
ENSG00000223972  0.000  0.117  0.000  0.000  0.000  0.058  0.000  0.000  0.000
ENSG00000227232  0.000  0.117  0.000  0.950  1.989  0.173  0.146  0.656  0.353
                HEA062 HEA063 HEA064 HEA065 HEA066 HEA067 HEA068 HEA069 HEA070
ENSG00000223972  0.000  0.000  0.062  0.000  0.000  0.000  0.000  0.000  0.000
ENSG00000227232  0.000  0.000  0.062  0.103  0.055  0.000  0.067  0.067  0.049
                HEA071 HEA181 HEA182 HEA183 HEA184 HEA185 HEA186
ENSG00000223972  0.000    0.0   0.00  0.000  0.000  0.000  0.000
ENSG00000227232  0.112    0.2   0.28  0.055  0.000  0.000  0.415
 [ reached getOption("max.print") -- omitted 4 rows ]
round(head(counts.keep), 3)
                COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155
ENSG00000227232  0.066  0.000  0.222  0.848  0.073  0.151  1.388  0.420  1.163
ENSG00000278267  0.199  0.072  0.148  0.000  0.220  0.000  0.000  0.070  0.174
                COV156 COV157 COV175 COV176 COV177 COV178 COV179 COV180 HEA061
ENSG00000227232  0.000  0.117  0.000  0.950  1.989  0.173  0.146  0.656  0.353
ENSG00000278267  0.128  0.117  0.224  0.219  0.435  0.000  0.219  0.000  0.235
                HEA062 HEA063 HEA064 HEA065 HEA066 HEA067 HEA068 HEA069 HEA070
ENSG00000227232  0.000  0.000  0.062  0.103  0.055  0.000  0.067  0.067  0.049
ENSG00000278267  0.288  0.193  0.495  0.154  0.493  0.427  0.134  0.469  0.098
                HEA071 HEA181 HEA182 HEA183 HEA184 HEA185 HEA186
ENSG00000227232  0.112  0.200   0.28  0.055  0.000  0.000  0.415
ENSG00000278267  0.337  0.200   0.14  0.328  0.699  0.304  0.178
 [ reached getOption("max.print") -- omitted 4 rows ]

3.3 Uso de clases específicas para manejar los datos

Cuando se trabaja con distintos objetos referidos a unos mismos datos, como la matriz de contajes y el objeto “targets”, es útil disponer de contenedores que permitan trabajar con todos ellos a la vez, lo que no sólo facilita el trabajo sino que ayuda a evitar “desincronizaciones”.

Éste es el caso de la clase ExpressionSet habitualmente utilizada con microarrays o de la clase que la generaliza, llamada SummarizedExperiment.

Para datos de contaje es habitual usar una clase similar a ExpressionSet llamada DGEList” pensadas para manejar datos de contajes , definida en el paquete edgeR. Esta clase, más simple que las anteriores, utiliza listas para almacenar recuentos de “reads” e información asociada de tecnologías de secuenciación. Puede encontrarse información al respecto en la ayuda del paquete edgeR.

dgeObj <- DGEList(counts = counts.keep, lib.size = colSums(counts.keep), norm.factors = rep(1,
    ncol(counts.keep)), samples = selectedTargets, group = selectedTargets$group,
    genes = rownames(counts.keep), remove.zeros = FALSE)
show(dgeObj)
An object of class "DGEList"
$counts
                    COV145     COV147    COV149     COV150     COV151
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991
ENSG00000278267 0.19934902 0.07170342 0.1479288 0.00000000 0.21962973
                    COV152    COV153     COV154   COV155    COV156     COV157
ENSG00000227232 0.15060933 1.3880853 0.41957493 1.163313 0.0000000 0.11736232
ENSG00000278267 0.00000000 0.0000000 0.06992916 0.174497 0.1279917 0.11736232
                   COV175    COV176    COV177    COV178     COV179    COV180
ENSG00000227232 0.0000000 0.9503052 1.9885554 0.1726206 0.14628246 0.6558384
ENSG00000278267 0.2236183 0.2193012 0.4349965 0.0000000 0.21942369 0.0000000
                   HEA061    HEA062    HEA063     HEA064     HEA065     HEA066
ENSG00000227232 0.3528731 0.0000000 0.0000000 0.06183451 0.10298192 0.05475139
ENSG00000278267 0.2352487 0.2880235 0.1926772 0.49467608 0.15447288 0.49276252
                   HEA067     HEA068     HEA069     HEA070    HEA071    HEA181
ENSG00000227232 0.0000000 0.06713275 0.06693886 0.04887949 0.1122282 0.2003508
ENSG00000278267 0.4268581 0.13426549 0.46857204 0.09775898 0.3366846 0.2003508
                    HEA182     HEA183    HEA184    HEA185     HEA186
ENSG00000227232 0.28009698 0.05465108 0.0000000 0.0000000 0.41517606
ENSG00000278267 0.14004849 0.32790646 0.6985910 0.3036215 0.17793260
 [ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...

$samples
       group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8            1  COV145   COVID  red
COV147 COVID 999530.1            1  COV147   COVID  red
COV149 COVID 999667.8            1  COV149   COVID  red
COV150 COVID 999872.5            1  COV150   COVID  red
COV151 COVID 999904.2            1  COV151   COVID  red
29 more rows ...

$genes
                          genes
ENSG00000227232 ENSG00000227232
ENSG00000278267 ENSG00000278267
ENSG00000233750 ENSG00000233750
ENSG00000268903 ENSG00000268903
ENSG00000269981 ENSG00000269981
33034 more rows ...

Uno de los aspectos interesantes de estas clases es la posibilidad de extraer partes de todos los objetos a la vez con el operador de “subsetting”.

dim(dgeObj)
[1] 33039    34
dgeObjShort <- dgeObj[, c(1:5, 11:15)]
# Library size information is stored in the samples slot
dgeObjShort$samples
       group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8            1  COV145   COVID  red
COV147 COVID 999530.1            1  COV147   COVID  red
COV149 COVID 999667.8            1  COV149   COVID  red
COV150 COVID 999872.5            1  COV150   COVID  red
COV151 COVID 999904.2            1  COV151   COVID  red
COV157 COVID 999582.4            1  COV157   COVID  red
COV175 COVID 998800.7            1  COV175   COVID  red
COV176 COVID 999667.6            1  COV176   COVID  red
COV177 COVID 999807.4            1  COV177   COVID  red
COV178 COVID 999721.9            1  COV178   COVID  red
colnames(dgeObjShort$counts)
 [1] "COV145" "COV147" "COV149" "COV150" "COV151" "COV157" "COV175" "COV176"
 [9] "COV177" "COV178"
# https://hbctraining.github.io/DGE_workshop_salmon/lessons/AnnotationDbi_lesson.html
annotations_orgDb <- AnnotationDbi::select(org.Hs.eg.db, # database
                                           keys = dgeObj$genes[,1],  # data to use for retrieval
                                           columns = c("SYMBOL", "ENTREZID","GENENAME"), # information to retreive for given data
                                           keytype = "ENSEMBL")

length(which(is.na(annotations_orgDb$SYMBOL) == FALSE))
[1] 23304
non_duplicates_idx <- which(duplicated(annotations_orgDb$SYMBOL == FALSE))
annotations_orgDb<- annotations_orgDb[non_duplicates_idx, ]
# https://hbctraining.github.io/Intro-to-R-flipped/lessons/09_reordering-to-match-datasets.html
cbind(dgeObj$genes[, 1], annotations_orgDb$ENSEMBL)[1:100, ]
       [,1]              [,2]             
  [1,] "ENSG00000227232" "ENSG00000278267"
  [2,] "ENSG00000278267" "ENSG00000233750"
  [3,] "ENSG00000233750" "ENSG00000269981"
  [4,] "ENSG00000268903" "ENSG00000241860"
  [5,] "ENSG00000269981" "ENSG00000279928"
  [6,] "ENSG00000241860" "ENSG00000279457"
  [7,] "ENSG00000279928" "ENSG00000228463"
  [8,] "ENSG00000279457" "ENSG00000237094"
  [9,] "ENSG00000228463" "ENSG00000230021"
 [10,] "ENSG00000237094" "ENSG00000225972"
 [11,] "ENSG00000230021" "ENSG00000225630"
 [12,] "ENSG00000225972" "ENSG00000237973"
 [13,] "ENSG00000225630" "ENSG00000229344"
 [14,] "ENSG00000237973" "ENSG00000248527"
 [15,] "ENSG00000229344" "ENSG00000198744"
 [16,] "ENSG00000248527" "ENSG00000228327"
 [17,] "ENSG00000198744" "ENSG00000237491"
 [18,] "ENSG00000228327" "ENSG00000230092"
 [19,] "ENSG00000237491" "ENSG00000177757"
 [20,] "ENSG00000230092" "ENSG00000228794"
 [21,] "ENSG00000177757" "ENSG00000225880"
 [22,] "ENSG00000228794" "ENSG00000288531"
 [23,] "ENSG00000225880" "ENSG00000234711"
 [24,] "ENSG00000288531" "ENSG00000272438"
 [25,] "ENSG00000234711" "ENSG00000230699"
 [26,] "ENSG00000272438" "ENSG00000223764"
 [27,] "ENSG00000230699" "ENSG00000187634"
 [28,] "ENSG00000223764" "ENSG00000188976"
 [29,] "ENSG00000187634" "ENSG00000187961"
 [30,] "ENSG00000188976" "ENSG00000187583"
 [31,] "ENSG00000187961" "ENSG00000187642"
 [32,] "ENSG00000187583" "ENSG00000272512"
 [33,] "ENSG00000187642" "ENSG00000188290"
 [34,] "ENSG00000272512" "ENSG00000187608"
 [35,] "ENSG00000188290" "ENSG00000224969"
 [36,] "ENSG00000187608" "ENSG00000188157"
 [37,] "ENSG00000224969" "ENSG00000217801"
 [ reached getOption("max.print") -- omitted 63 rows ]
reorder_idx <- match(dgeObj$genes[, 1], annotations_orgDb$ENSEMBL)
annotations_orgDb_reordered <- annotations_orgDb[reorder_idx, ]
cbind(dgeObj$genes[, 1], annotations_orgDb_reordered$ENSEMBL)[1:100, ]
       [,1]              [,2]             
  [1,] "ENSG00000227232" NA               
  [2,] "ENSG00000278267" "ENSG00000278267"
  [3,] "ENSG00000233750" "ENSG00000233750"
  [4,] "ENSG00000268903" NA               
  [5,] "ENSG00000269981" "ENSG00000269981"
  [6,] "ENSG00000241860" "ENSG00000241860"
  [7,] "ENSG00000279928" "ENSG00000279928"
  [8,] "ENSG00000279457" "ENSG00000279457"
  [9,] "ENSG00000228463" "ENSG00000228463"
 [10,] "ENSG00000237094" "ENSG00000237094"
 [11,] "ENSG00000230021" "ENSG00000230021"
 [12,] "ENSG00000225972" "ENSG00000225972"
 [13,] "ENSG00000225630" "ENSG00000225630"
 [14,] "ENSG00000237973" "ENSG00000237973"
 [15,] "ENSG00000229344" "ENSG00000229344"
 [16,] "ENSG00000248527" "ENSG00000248527"
 [17,] "ENSG00000198744" "ENSG00000198744"
 [18,] "ENSG00000228327" "ENSG00000228327"
 [19,] "ENSG00000237491" "ENSG00000237491"
 [20,] "ENSG00000230092" "ENSG00000230092"
 [21,] "ENSG00000177757" "ENSG00000177757"
 [22,] "ENSG00000228794" "ENSG00000228794"
 [23,] "ENSG00000225880" "ENSG00000225880"
 [24,] "ENSG00000288531" "ENSG00000288531"
 [25,] "ENSG00000234711" "ENSG00000234711"
 [26,] "ENSG00000272438" "ENSG00000272438"
 [27,] "ENSG00000230699" "ENSG00000230699"
 [28,] "ENSG00000223764" "ENSG00000223764"
 [29,] "ENSG00000187634" "ENSG00000187634"
 [30,] "ENSG00000188976" "ENSG00000188976"
 [31,] "ENSG00000187961" "ENSG00000187961"
 [32,] "ENSG00000187583" "ENSG00000187583"
 [33,] "ENSG00000187642" "ENSG00000187642"
 [34,] "ENSG00000272512" "ENSG00000272512"
 [35,] "ENSG00000188290" "ENSG00000188290"
 [36,] "ENSG00000187608" "ENSG00000187608"
 [37,] "ENSG00000224969" "ENSG00000224969"
 [ reached getOption("max.print") -- omitted 63 rows ]
dgeObj$genes <- annotations_orgDb_reordered
dgeObj$genes[1:100, ]
             ENSEMBL       SYMBOL  ENTREZID
NA              <NA>         <NA>      <NA>
2    ENSG00000278267    MIR6859-1 102466751
3    ENSG00000233750       CICP27 100420257
NA.1            <NA>         <NA>      <NA>
5    ENSG00000269981         <NA>      <NA>
6    ENSG00000241860         <NA>      <NA>
7    ENSG00000279928         <NA>      <NA>
8    ENSG00000279457       WASH9P 102723897
9    ENSG00000228463    RPL23AP21    728481
10   ENSG00000237094         <NA>      <NA>
11   ENSG00000230021 LOC101928626 101928626
12   ENSG00000225972     MTND1P23 100887749
13   ENSG00000225630     MTND2P28 100652939
14   ENSG00000237973     MTCO1P12 107075141
15   ENSG00000229344     MTCO2P12 107075310
16   ENSG00000248527     MTATP6P1 106480796
17   ENSG00000198744     MTCO3P12 107075270
18   ENSG00000228327         <NA>      <NA>
                                            GENENAME
NA                                              <NA>
2                                    microRNA 6859-1
3    capicua transcriptional repressor pseudogene 27
NA.1                                            <NA>
5                                               <NA>
6                                               <NA>
7                                               <NA>
8           WAS protein family homolog 9, pseudogene
9               ribosomal protein L23a pseudogene 21
10                                              <NA>
11                      uncharacterized LOC101928626
12                              MT-ND1 pseudogene 23
13                              MT-ND2 pseudogene 28
14                              MT-CO1 pseudogene 12
15                              MT-CO2 pseudogene 12
16                              MT-ATP6 pseudogene 1
17                              MT-CO3 pseudogene 12
18                                              <NA>
 [ reached 'max' / getOption("max.print") -- omitted 82 rows ]

Aunque podríamos haber creado el objeto a partir de todas las muestras, y haber realizado la extracción de genes y muestras posteriormente, hemos optado por no hacerlo para facilitar el seguimiento del proceso.

3.4 Normalización

Además de estandarizar los contajes, es importante eliminar otros sesgos de composición entre librerías. Esto puede hacerse aplicando la normalización por el método TMM que genera un conjunto de factores de normalización, tal que producto de estos factores y los tamaños de librería (el número de secuencias de cada muestra) definen el tamaño efectivo de dichas muestras, es decir el peso real que se les asignará en las comparaciones posteriores.

Aunque esto puede parecer artificial, no lo es porque la normalización tiene en cuenta otros factores, como el sesgo de composición entre librerías, que podrían hacer que los mismos valores en distintas muestras no reflejaran su importancia relativa.

La función calcNormFactors, de la librería edgeR, calcula los factores de normalización mencionados.

library(edgeR)
dgeObj_norm <- calcNormFactors(dgeObj)
head(dgeObj_norm$samples, 10)
       group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8    1.0018278  COV145   COVID  red
COV147 COVID 999530.1    1.1010464  COV147   COVID  red
COV149 COVID 999667.8    0.9553476  COV149   COVID  red
COV150 COVID 999872.5    0.9479175  COV150   COVID  red
COV151 COVID 999904.2    0.8996072  COV151   COVID  red
COV152 COVID 999763.7    1.0033710  COV152   COVID  red
COV153 COVID 999874.5    0.9870834  COV153   COVID  red
COV154 COVID 999736.9    0.8538699  COV154   COVID  red
COV155 COVID 999399.0    0.8230210  COV155   COVID  red
COV156 COVID 999868.7    0.8878678  COV156   COVID  red

Esto no modificará la matriz de contajes, pero actualizará los factores de normalización en el objeto DGEList (sus valores predeterminados son 1).

head(dgeObj$samples)
       group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8            1  COV145   COVID  red
COV147 COVID 999530.1            1  COV147   COVID  red
COV149 COVID 999667.8            1  COV149   COVID  red
COV150 COVID 999872.5            1  COV150   COVID  red
COV151 COVID 999904.2            1  COV151   COVID  red
COV152 COVID 999763.7            1  COV152   COVID  red
head(dgeObj_norm$samples)
       group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8    1.0018278  COV145   COVID  red
COV147 COVID 999530.1    1.1010464  COV147   COVID  red
COV149 COVID 999667.8    0.9553476  COV149   COVID  red
COV150 COVID 999872.5    0.9479175  COV150   COVID  red
COV151 COVID 999904.2    0.8996072  COV151   COVID  red
COV152 COVID 999763.7    1.0033710  COV152   COVID  red

Es decir, aunque no se observen cambios en la matriz de contajes, cuando se utilizan estos factores de normalización en algún cálculo la importancia de las distintas columnas se tendrá en cuenta.

3.5 Transformación logarítmica

Las transformaciones anteriores buscan compensar el tamaño distinto de las librerías o la distinta composición de éstas, pero las distribuciones de los contajes en cada muestra son asimétricas.

boxplot(dgeObj_norm$counts, col = dgeObj_norm$samples$cols, las = 2, cex.axis = 0.7,
    main = "Contajes normalizados", ylim = c(0, 10000))

Para finalizar el preprocesado se toman logaritmo base dos de los contajes.

log2count_norm <- cpm(dgeObj_norm, log = TRUE)
source("./Rcode/niceBoxPlot.R")
gg1 <- niceBoxPlot(log2count_norm)
ggsave(filename = "figures/boxplot.png", gg1, dpi = 300, bg = "transparent")
gg1

También se puede hacer usando plotly.

niceBoxPlotly(log2count_norm)

Esta será nuestra matriz de partida para los análisis siguientes

4 Exploración de los datos

Una vez descartados los genes poco expresados y con los recuentos almacenados en un objeto DGEList, podemos`proceder a realizar algunos gráficos exploratorios para determinar si los datos aparentan buena calidad y/o si presentan algun problema.

4.1 Distribución de los contajes

Un diagrama de cajas con los datos, normalizados o no, muestra que la distribución de los contajes es muy asimétrica, lo que justifica la decisión de trabajar con los logaritmos de los datos.

La transformación logarítmica puede hacerse directamente pero es mejor usar la función cpm, como se ha hecho, que agrega una pequeña cantidad para evitar tomar logaritmos de cero.

par(mfrow = c(2, 1))
rawCounts <- dgeObj_norm$counts
boxplot(rawCounts, ylab = "CPM", las = 2, xlab = "", col = dgeObj$samples$cols, cex.axis = 0.6,
    main = "Distribución de contajes")
boxplot(log2count_norm, ylab = "Log2-CPM", las = 2, xlab = "", col = dgeObj$samples$cols,
    cex.axis = 0.6, main = "Distribución de log(contajes)")
abline(h = median(log2count_norm), col = "blue")

4.2 Análisis de similaridad entre las muestras

4.2.1 Distancia entre muestras

La función dist permite calcular una matriz de distancias que contiene las comparaciones dos a dos entre todas las muestras. Por defecto se utiliza una distancia euclídea.

sampleDists <- dist(t(log2count_norm))
round(sampleDists, 1)
       COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155 COV156
COV147   80.0                                                               
COV149   70.7   82.0                                                        
       COV157 COV175 COV176 COV177 COV178 COV179 COV180 HEA061 HEA062 HEA063
COV147                                                                      
COV149                                                                      
       HEA064 HEA065 HEA066 HEA067 HEA068 HEA069 HEA070 HEA071 HEA181 HEA182
COV147                                                                      
COV149                                                                      
       HEA183 HEA184 HEA185
COV147                     
COV149                     
 [ reached getOption("max.print") -- omitted 31 rows ]

Las matrices de distancias se pueden visualizar directamente con un heatmap.

# par(mfrow=c(1,1))
library(factoextra)
fviz_dist(sampleDists)

Como puede verse _las muestras tienden a agruparse por el factor SANO/COVID, aunque una de las muestras. COV155 se separa del resto de las del grupo COVID.

4.2.2 Agrupamiento jerárquico

Un agrupamiento jerárquico proporciona una representación alternativa, también basada en la matriz de distancias.

hc <- hclust(sampleDists)
plot(hc, labels = colnames(log2count_norm), main = "Agrpamiento jerárquico de las muestras",
    cex = 0.8)

El siguiente es un gráfico mejor.

source("./Rcode/niceDendrogram.R")
library("ggdendro", "dendextend")
hc <- hclust(sampleDists, method = "ward.D2")
gg2 <- niceDendrogram(hc)
ggsave("figures/dendograma.png", gg2, bg = "transparent")  # , dpi = 300
gg2

Una de las muestras COVID parece más similar a las saludables que a las otras COVID.

4.2.3 Visualización en dimensión reducida

col.status <- dgeObj_norm$samples$cols
plotMDS(log2count_norm, col = col.status, main = "Status", cex = 0.7)

Con el paquete emojifont se pueden usar emoticonos para distinguir los pacientes sanos de los que tienen COVID.

source("./Rcode/niceMDS.R")
gg3 <- niceMDS(dgeObj_norm)
ggplot2::ggsave("figures/MDS_emotis.png", gg3, bg = "transparent")  # ,dpi=300
gg3

Como puede verse, el gráfico muestra la misma agrupación “natural” y el mismo comportamiento atípico de una muestra.

5 Análisis de expresión diferencial

El objetivo del análisis de expresión diferencial es seleccionar genes cuya expresión difiere entre grupos.

5.1 Selección de genes usando limma-Voom

La ventaja principal de esta aproximación es que permite trabajar con toda la flexibilidad de los modelos lineales para representar diseños experimentales, y, en muchos casos , aprovechar la experiencia previa del usuario en el manejo de limma.

5.1.1 Matriz de diseño y de contrastes

Utilizando la variable group podemos definir una matriz de diseño y, sobre ésta, los contrastes que nos interesan.

group = as.factor(dgeObj_norm$samples$group)
design = model.matrix(~0 + group)
colnames(design) = gsub("group", "", colnames(design))
row.names(design) = sampleNames
design
       COVID SANO
COV145     1    0
COV147     1    0
COV149     1    0
COV150     1    0
COV151     1    0
COV152     1    0
COV153     1    0
COV154     1    0
COV155     1    0
COV156     1    0
COV157     1    0
COV175     1    0
COV176     1    0
COV177     1    0
COV178     1    0
COV179     1    0
COV180     1    0
HEA061     0    1
HEA062     0    1
HEA063     0    1
HEA064     0    1
HEA065     0    1
HEA066     0    1
HEA067     0    1
HEA068     0    1
HEA069     0    1
HEA070     0    1
HEA071     0    1
HEA181     0    1
HEA182     0    1
HEA183     0    1
HEA184     0    1
HEA185     0    1
HEA186     0    1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

Dado que estamos interesados en las diferencias entre los grupos, necesitamos especificar qué comparaciones queremos llevar a cabo. Las comparaciones de interés se puede especificar utilizando la función makeContrasts. La matriz de contraste indica qué columnas de la matriz design vamos a comparar. En este caso tan sólo se llevará a cabo una comparación.

cont.matrix = makeContrasts(CONTROLvsCOVID = COVID - SANO, levels = colnames(design))
cont.matrix
       Contrasts
Levels  CONTROLvsCOVID
  COVID              1
  SANO              -1

5.1.2 Transformación de los contajes

voomObj <- voom(dgeObj_norm, design)
voomObj
An object of class "EList"
$genes
                          genes
ENSG00000227232 ENSG00000227232
ENSG00000278267 ENSG00000278267
ENSG00000233750 ENSG00000233750
ENSG00000268903 ENSG00000268903
ENSG00000269981 ENSG00000269981
33034 more rows ...

$targets
       group  lib.size norm.factors samples group.1 cols
COV145 COVID 1001653.3    1.0018278  COV145   COVID  red
COV147 COVID 1100529.0    1.1010464  COV147   COVID  red
COV149 COVID  955030.2    0.9553476  COV149   COVID  red
COV150 COVID  947796.6    0.9479175  COV150   COVID  red
COV151 COVID  899521.0    0.8996072  COV151   COVID  red
29 more rows ...

$E
                    COV145     COV147     COV149      COV150     COV151
ENSG00000227232 -0.8223650 -1.1381984 -0.4037625  0.50837551 -0.6500950
ENSG00000278267 -0.5183002 -0.9448596 -0.5597126 -0.92265092 -0.3219038
                    COV152     COV153     COV154     COV155     COV156
ENSG00000227232 -0.6246522  0.9358596  0.1073289  1.0159241 -0.8282289
ENSG00000278267 -1.0045156 -0.9810644 -0.5828556 -0.2862517 -0.4994115
                     COV157     COV175     COV176     COV177     COV178
ENSG00000227232 -0.74055345 -0.9859431  0.5497283  1.3062240 -0.6882933
ENSG00000278267 -0.74055345 -0.4526424 -0.4619603 -0.1060517 -1.1161581
                      COV179     COV180      HEA061      HEA062     HEA063
ENSG00000227232 -0.654532686  0.1342634 -0.20656331 -1.01143480 -1.1012204
ENSG00000278267 -0.499855856 -1.0746764 -0.42066201 -0.35512422 -0.6309653
                    HEA064     HEA065      HEA066     HEA067     HEA068
ENSG00000227232 -0.9265423 -0.7259291 -0.78574256 -1.1218171 -0.9590727
ENSG00000278267 -0.1024608 -0.6077105  0.05386471 -0.2313968 -0.7976723
                    HEA069     HEA070     HEA071     HEA181     HEA182
ENSG00000227232 -1.0099326 -0.8215555 -0.6475041 -0.4772445 -0.4479095
ENSG00000278267 -0.2372664 -0.6984811 -0.1968898 -0.4772445 -0.7333818
                    HEA183     HEA184     HEA185     HEA186
ENSG00000227232 -0.8517698 -1.0495769 -0.9932084 -0.1628122
ENSG00000278267 -0.2738825  0.2117625 -0.3086203 -0.5957197
 [ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...

$weights
          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
[1,]  6.976324  6.423290  7.345017  7.414329  7.913732  6.966953  7.070951
[2,] 13.344117  9.792171 16.233083 16.785098 21.291716 13.276959 14.091120
          [,8]      [,9]     [,10]     [,11]     [,12]     [,13]     [,14]
[1,]  8.568175  9.129126  8.046130  6.792973  7.049157  7.046103  6.946915
[2,] 27.856177 34.175509 22.694541 12.084889 13.891449 13.863695 13.134382
         [,15]     [,16]     [,17]     [,18]     [,19]     [,20]     [,21]
[1,]  6.497274  6.878669  6.667025 21.376991 19.139568 14.546448 14.827401
[2,] 10.280091 12.659254 11.283734  8.528972  8.192672  7.560903  7.603134
         [,22]     [,23]     [,24]     [,25]     [,26]     [,27]     [,28]
[1,] 20.095913 24.685105 13.715056 13.110198 11.661751 22.943084 24.332783
[2,]  8.340050  8.955089  7.428218  7.308304  7.031769  8.740787  8.912791
         [,29]     [,30]     [,31]     [,32]     [,33]     [,34]
[1,] 22.360747 15.055081 19.758022 16.994284 20.284112 17.780471
[2,]  8.666115  7.636840  8.288627  7.906856  8.368394  8.008534
 [ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...

$design
       COVID SANO
COV145     1    0
COV147     1    0
COV149     1    0
COV150     1    0
COV151     1    0
29 more rows ...

5.1.3 Selección de genes diferencialmente expresados

Como en el caso de los microarrays el objeto voomObj y las matrices de diseño y contrastes se utilizaran para ajustar un modelo y, a continuación realizar las comparaciones especificadas sobre el modelo ajustado. El proceso finaliza con la regularización del estimador del error usando la función eBayes.

fit <- lmFit(voomObj)
fit.cont <- contrasts.fit(fit, cont.matrix)
fit.cont <- eBayes(fit.cont)

5.1.4 Top tables

Los resultados de un análisis de expresión diferencial se pueden extraer con la función topTable. Esta función genera una tabla de resultados cuyas columnas contienen información acerca de los genes y la diferencia entre los grupos comparados. Concretamente:

toptab <- topTable(fit.cont, coef = 1, sort.by = "p", number = nrow(fit.cont))
head(toptab)
                          genes    logFC  AveExpr        t      P.Value
ENSG00000167900 ENSG00000167900 2.832423 4.111065 12.36669 8.088344e-15
ENSG00000090889 ENSG00000090889 2.476447 1.292066 11.83766 3.015984e-14
ENSG00000111206 ENSG00000111206 2.423923 2.018435 11.80083 3.309678e-14
ENSG00000166851 ENSG00000166851 3.223189 2.615540 11.65570 4.781203e-14
ENSG00000135476 ENSG00000135476 2.236975 1.705254 11.53985 6.425179e-14
ENSG00000123485 ENSG00000123485 2.553959 1.914888 11.46276 7.828950e-14
                   adj.P.Val        B
ENSG00000167900 2.672308e-10 23.54910
ENSG00000090889 3.644948e-10 22.09352
ENSG00000111206 3.644948e-10 21.99541
ENSG00000166851 3.949154e-10 21.70990
ENSG00000135476 4.013294e-10 21.35292
ENSG00000123485 4.013294e-10 21.18953

5.1.5 Visualización de los resultados

5.1.5.1 Significación biológica vs. Significación Estadística

Para visualizar los resultados podemos usar un volcanoPlot:

volcanoplot(fit.cont, coef = 1, highlight = 100, main = "COVID vs SANO")

Es posible también crear Volcanoplots interactivos

source("./Rcode/niceInteractiveVolcano.R")
niceInteractiveVolcano(fit.cont)

O volcano plots más sofisticados

library("gridExtra")
library("plotly")
diseased_vs_healthy <- data.frame(fit.cont$genes[, 1], fit.cont$coefficients[, 1],
    fit.cont$p.value[, 1])
colnames(diseased_vs_healthy) <- c("geneid", "foldchange", "adjpval")
diseased_vs_healthy <- diseased_vs_healthy %>%
    mutate(gene_type = case_when(foldchange >= 2 & adjpval <= 0.05 ~ "up", foldchange <=
        -2 & adjpval <= 0.05 ~ "down", TRUE ~ "ns"))

# Obtain gene_type counts
# ------------------------------------------------------
diseased_vs_healthy %>%
    count(gene_type)
  gene_type     n
1      down     4
2        ns 32837
3        up   198
p<-  ggplot(diseased_vs_healthy,aes(x = foldchange,
             y = -log10(adjpval),
             fill = gene_type,    
             size = gene_type,
             alpha = gene_type)) + 
  geom_point(shape = 21, # Specify shape and colour as fixed local parameters    
             colour = "black") + 
  geom_hline(yintercept = -log10(0.05),
             linetype = "dashed") + 
  geom_vline(xintercept = c(-2, 2),
             linetype = "dashed") +
  scale_fill_manual(values = cols) + # Modify point colour
  scale_size_manual(values = sizes) + # Modify point size
  scale_alpha_manual(values = alphas) + # Modify point transparency
  scale_x_continuous(breaks = c(seq(-6.5, 6.5, 1)),       
                     limits = c(-7, 7))  
p1<-p+theme(
         panel.background = element_rect(fill='transparent'),
         plot.background = element_rect(fill='transparent', color=NA),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.key = element_rect(colour = NA, fill = NA),
         legend.background = element_rect(fill='transparent'),
         legend.box.background = element_rect(fill='transparent'))+
  annotate("text", x=6.9, y=7.8, label= "IFI27")+
  annotate("text", x=5.8, y=11.7, label= "IGHG1") +
  annotate("text", x=5.7, y=10.5, label= "IGHG1-24") + 
  annotate("text", x=5.5, y=4.3, label= "ALAS2") + 
  annotate("text", x=5.2, y=2.9, label= "HBB") + 
  annotate("text", x=-1.0, y=11.4, label= "CACNA2D3") + 
  annotate("text", x=-1.5, y=7.4, label= "UICLM") + 
  annotate("text", x=-1.45, y=5.2, label= "HLA-DQB2") + 
  annotate("text", x=-1.4, y= 4.0, label= "HLA-DQA1") 


#display boxplot
p1

ggsave("figures/volcanoplot1.png", p1, dpi = 300, bg = "transparent")
# add a grouping column; default value is 'not significant'

diff_df <- diseased_vs_healthy
colnames(diff_df) <- c("geneid", "foldchange", "adjpval", "genetype")
diff_df["group"] <- "NotSignificant"

# for our plot, we want to highlight FDR < 0.05 (significance level) Fold
# Change > 1.5

# change the grouping for the entries with significance but not a large enough
# Fold change
diff_df[which(diff_df["adjpval"] < 0.05 & abs(diff_df["foldchange"]) < 1.5), "group"] <- "Significant"

# change the grouping for the entries a large enough Fold change but not a low
# enough p value
diff_df[which(diff_df["adjpval"] > 0.05 & abs(diff_df["foldchange"]) > 1.5), "group"] <- "FoldChange"

# change the grouping for the entries with both significance and large enough
# fold change
diff_df[which(diff_df["adjpval"] < 0.05 & abs(diff_df["foldchange"]) > 1.5), "group"] <- "Significant&FoldChange"


# Find and label the top peaks..
top_peaks <- diff_df[with(diff_df, order(foldchange, adjpval)), ][1:5, ]
top_peaks <- rbind(top_peaks, diff_df[with(diff_df, order(-foldchange, adjpval)),
    ][1:5, ])

# Add gene labels for all of the top genes we found here we are creating an
# empty list, and filling it with entries for each row in the dataframe each
# list entry is another list with named items that will be used by Plot.ly
a <- list()
for (i in seq_len(nrow(top_peaks))) {
    m <- top_peaks[i, ]
    a[[i]] <- list(x = m[["foldchange"]], y = -log10(m[["adjpval"]]), text = m[["geneid"]],
        xref = "x", yref = "y", showarrow = TRUE, arrowhead = 0.5, ax = 20, ay = -40)
}
# make the Plot.ly plot
p <- plot_ly(data = diff_df, x = ~foldchange, y = ~-log10(adjpval), text = ~geneid,
    mode = "markers", color = ~group, xaxis = ~c(-7, 7)) %>%
    layout(title = "") %>%
    layout(annotations = a) %>%
    layout(xaxis = list(range = c(-7, 7)))

p

5.1.5.2 Perfiles de expresión

Con el fin de observar si existen perfiles de expresión diferenciados podemo realizar un mapa de colores con los genes más diferencialmente expresados.

Es decir, fijamos un criterio de selección de genes y retenemos aquellos componentes de la tabla de resultados que lo cumplen. Por ejemplo: Genes con un p-balor ajustado inferior a 0.001 y un `fold-change’ superior a 2.

topGenesBas <- rownames(subset(toptab, (abs(logFC) > 2) & (adj.P.Val < 0.01)))
length(topGenesBas)
[1] 199

Con la matriz de expresión de los genes que verifican dicha condición se puede construir un heatmap.

library(pheatmap)
mat <- log2count_norm[topGenesBas, ]
mat <- mat - rowMeans(mat)
pheatmap(mat, fontsize_col = 14, fontsize_row = 1)

png("figures/pheatmap_transp.png", units = "in", width = 30, height = 22, res = 600,
    bg = "transparent")
par(bg = NA)
pheatmap(mat, fontsize = 30, fontsize_col = 30, fontsize_row = 1, treeheight_row = 60,
    treeheight_col = 180)
dev.off()
png 
  3

También es posible crear un Heatmap interactivo:

source("./Rcode/niceInteractiveHeatMap.R")
mat <- log2count_norm[topGenesBas, ]
mat <- mat - rowMeans(mat)
niceInteractiveHeatMap(mat)

Los dos grupos estan diferenciados, sobre todo en un subconjunto de genes en donde las expresiones toman signos (colores) distintos. En otro de los grupos parece que, en el grupo de individuos sanos los genes diferencialmente expresados se encuentran sobre-expresados en el grupo COVID, y apenas expresados en el grupo sanos.

5.2 Análisis de expresión diferencial usando el paquete edgeR

El análisis con edgeR es similar al anterior (se originan en el mismo equipo de investigación) pero la modelización es distinta.

El análisis utiliza un GLM pero, en una forma que recuerda a lo que se hace con limma, realiza un paso adicional,en el que se calcula una estimación mejorada de la dispersión (variabilidad) de las muestras que integra las estimaciones individuales y la global mediante estimación Bayes empírica.

y = estimateDisp(dgeObj_norm, design, robust = TRUE)
plotBCV(y)

Con este objeto, que añade a los contajes normalizados los estimadores mejorados de dispersión, se ajusta un modelo lineal generalizado con distribución binomial para los errores.

fit <- glmQLFit(y, design, robust = TRUE)

Una vez ajustado el modelo se procede a construir el contraste y realizar el test.

De hecho podemos usar la matriz de contrastes que construímos para limma voom, en la misma forma que hemos reutilizado la de diseño.

res <- glmQLFTest(fit, contrast = cont.matrix)

Los resultados se almacenan en un objeto similar a la ’ topTable’ de limma.

topTags_edge <- topTags(res, n = dim(log2count_norm)[1])  # todos los genes
head(topTags_edge)
Coefficient:  1*COVID -1*SANO 
                          genes    logFC   logCPM        F       PValue
ENSG00000123485 ENSG00000123485 2.969009 2.874508 143.3978 6.410064e-15
ENSG00000167900 ENSG00000167900 3.052165 4.946915 165.9368 1.024164e-14
ENSG00000187837 ENSG00000187837 1.719557 4.813853 135.1654 1.889383e-14
ENSG00000111206 ENSG00000111206 2.834520 2.919472 132.0747 2.767000e-14
ENSG00000178999 ENSG00000178999 2.987952 3.225967 135.0640 4.056862e-14
ENSG00000166851 ENSG00000166851 3.638997 3.745942 144.7154 5.207210e-14
                         FDR
ENSG00000123485 1.691868e-10
ENSG00000167900 1.691868e-10
ENSG00000187837 2.080777e-10
ENSG00000111206 2.285473e-10
ENSG00000178999 2.680693e-10
ENSG00000166851 2.867350e-10

Podemos seleccionar los genes más diferencialmente expresados de la misma forma que hicimos con limma-voom

topGenes_edge <- rownames(subset(topTags_edge$table, (abs(logFC) > 2) & (FDR < 0.01)))
length(topGenes_edge)
[1] 376

Obsérvese que la lista de genes seleccionados es muy similar en ambos casos.

library(ggvenn)
x = list(LimmaVoom = topGenesBas, edgeR = topGenes_edge)
ggvenn(x, fill_color = c("#0073C2FF", "#EFC000FF"), stroke_size = 0.5, set_name_size = 3)

6 Anotación de resultados y análisis de significación biológica

Para el análisis de significación se utilizan dos listas de transcritos:

topGenes <- union(topGenesBas, topGenes_edge)
length(topGenes)
[1] 383
universe <- rownames(toptab)
length(universe)
[1] 33039

6.1 Anotación de los identificadores

Esto es posible, y de hecho sencillo de llevar a cabo, usando el paquete annotate.

library(org.Hs.eg.db)
AnnotationDbi::keytypes(org.Hs.eg.db)
 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
 [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
[11] "GENETYPE"     "GO"           "GOALL"        "IPI"          "MAP"         
[16] "OMIM"         "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"        
[21] "PMID"         "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
[26] "UNIPROT"     
topAnots = AnnotationDbi::select(org.Hs.eg.db, topGenes, c("SYMBOL", "ENTREZID",
    "GENENAME"), keytype = "ENSEMBL")
head(topAnots)
          ENSEMBL SYMBOL ENTREZID                                   GENENAME
1 ENSG00000167900    TK1     7083                         thymidine kinase 1
2 ENSG00000090889  KIF4A    24137                   kinesin family member 4A
3 ENSG00000111206  FOXM1     2305                            forkhead box M1
4 ENSG00000166851   PLK1     5347                         polo like kinase 1
5 ENSG00000135476  ESPL1     9700 extra spindle pole bodies like 1, separase
6 ENSG00000123485  HJURP    55355      Holliday junction recognition protein
dim(topAnots)
[1] 385   4

Como puede verse, el número de anotaciones es el mismo que el de identificadores ENSEMBL, lo que podría llevar a pensar que, es posible que, antes de subir los datos a GEO se hayan agrupado los contajes por genes.

Para la anotación del universo se procederá igual.

univAnots = AnnotationDbi::select(org.Hs.eg.db, universe, c("SYMBOL", "ENTREZID",
    "GENENAME"), keytype = "ENSEMBL")
head(univAnots)
          ENSEMBL SYMBOL ENTREZID                                   GENENAME
1 ENSG00000167900    TK1     7083                         thymidine kinase 1
2 ENSG00000090889  KIF4A    24137                   kinesin family member 4A
3 ENSG00000111206  FOXM1     2305                            forkhead box M1
4 ENSG00000166851   PLK1     5347                         polo like kinase 1
5 ENSG00000135476  ESPL1     9700 extra spindle pole bodies like 1, separase
6 ENSG00000123485  HJURP    55355      Holliday junction recognition protein
dim(univAnots)
[1] 33204     4

En este caso se observa como hay más anotaciones que transcritos, lo que sugiere que múltiples transcritos han sido mapeados en el mismo gen.

6.2 Análisis de enriquecimiento

El paquete clusterProfiler admite identificadores de tipo ENSEMBL y permite gran variedad de análisis complementarios al enriquecimiento por lo que, es una de las mejores opciones para el análisis de significación biológica.

library(clusterProfiler)
library(org.Hs.eg.db)
ego = enrichGO(gene = topGenes, universe = universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db,
    ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05, readable = TRUE)
head(ego[, -8], n = 5)
                   ID
GO:0006958 GO:0006958
GO:0006910 GO:0006910
GO:0002455 GO:0002455
GO:0006956 GO:0006956
GO:0002377 GO:0002377
                                                              Description
GO:0006958                       complement activation, classical pathway
GO:0006910                                      phagocytosis, recognition
GO:0002455 humoral immune response mediated by circulating immunoglobulin
GO:0006956                                          complement activation
GO:0002377                                      immunoglobulin production
           GeneRatio   BgRatio       pvalue     p.adjust       qvalue Count
GO:0006958    52/340  92/16054 1.454551e-63 4.321470e-60 4.020685e-60    52
GO:0006910    50/340  86/16054 5.387595e-62 8.003272e-59 7.446223e-59    50
GO:0002455    53/340 105/16054 3.515964e-61 3.481976e-58 3.239622e-58    53
GO:0006956    52/340 110/16054 4.674208e-58 3.471768e-55 3.230124e-55    52
GO:0002377    63/340 198/16054 3.010824e-57 1.789032e-54 1.664510e-54    63

Con los resultados del análisis de enriquecimiento se pueden llevar a cabo distintas visualizaciones cuya interpretación exacta puede verse en el manual de clusterProfiler

dotplot(ego, showCategory = 7)

Un dotplot un poco más trabajado

## sense transparencia
png("figures/dotplot.png", units="in", width=14, height=9, res=600)
dotplot(ego, showCategory=15, font.size = 12)
invisible(dev.off())
## amb transparencia
ego2 = clusterProfiler::simplify(ego, cutoff = 0.5, by = "p.adjust")
p <- dotplot(ego2, showCategory=15, font.size = 20) + 
  scale_colour_gradient2(low = "red", mid = "orange", high = "orange", midpoint = 0.0010) +
  theme(panel.background  = element_rect(fill = "transparent"),             # bg of the panel
        plot.background   = element_rect(fill = "transparent", color = NA), # bg of the plot
        legend.background = element_rect(fill = "transparent"),             # get rid of legend bg
        legend.box.background = element_rect(fill = "transparent"))         # bg legend boxes
ggsave("figures/dotplot_transp.png", width=14, height=9, dpi=300, bg = "transparent")
p

library(clusterProfiler)
library(ggplot2)
ego2 = clusterProfiler::simplify(ego, cutoff = 0.5, by = "p.adjust")
png("figures/cnetplot_transp.png", units = "in", width = 24, height = 16, res = 600,
    bg = "transparent")
par(bg = NA)
a <- cnetplot(ego2, showCategory = 5, cex_category = 1, cex_label_category = 2.5,
    cex_gene = 1, cex_label_gene = 1, circular = FALSE, colorEdge = TRUE)
a
invisible(dev.off())
a

library(enrichplot)
goplot(ego2, showCategory = 6, cex = 0.1)

heatplot(ego2)

term_similarity_matrix = pairwise_termsim(ego)
emapplot(term_similarity_matrix, showCategory = 15, group_category = TRUE, group_legend = TRUE)

7 Referencias cortas

Referencias largas

Allison, David B., Xiangqin Cui, Grier P. Page, and Mahyar Sabripour. 2006. “Microarray data analysis: From disarray to consolidation and consensus.” Nature Reviews Genetics 7 (1): 55–65. https://doi.org/10.1038/nrg1749.
Arunachalam, Prabhu S et al. 2020. “Systems Biological Assessment of Immunity to Mild Versus Severe COVID-19 Infection in Humans.” Science 369 (6508): 1210–20. https://doi.org/10.1126/science.abc6261.
Benjamini, Yoav, and Yosef Hochberg. 1995. “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Source Journal of the Royal Statistical Society. Series B (Methodological) 57 (1): 289–300. http://www.jstor.org/stable/2346101 http://www.jstor.org/ http://www.jstor.org/action/showPublisher?publisherCode=black.
Bushel, Pierre. 2013. Pvca: Principal Variance Component Analysis (PVCA).
Carlson, Marc. 2017. Org.mm.eg.db: Genome Wide Annotation for Mouse.
Carvalho, Benilton. 2015. Pd.mogene.2.1.st: Platform Design Info for Affymetrix MoGene-2.1-St.
Carvalho, Benilton S, and Rafael A Irizarry. 2010. “A Framework for Oligonucleotide Microarray Preprocessing.” Bioinformatics 26 (19): 2363–67. https://doi.org/10.1093/bioinformatics/btq431.
Chrominski, Kornel, and Magdalena Tkacz. 2015. “Comparison of High-Level Microarray Analysis Methods in the Context of Result Consistency.” PLOS ONE 10 (6): e0128845. https://doi.org/10.1371/JOURNAL.PONE.0128845.
Dahl, David B. 2016. Xtable: Export Tables to LaTeX or HTML. https://CRAN.R-project.org/package=xtable.
Davis, Sean, and Paul Meltzer. 2007. “GEOquery: A Bridge Between the Gene Expression Omnibus (GEO) and BioConductor.” Bioinformatics 14: 1846–47.
Draghici, Sorin. 2012. Statistics and data analysis for microarrays using R and Bioconductor. CRC Press. https://www.crcpress.com/Statistics-and-Data-Analysis-for-Microarrays-Using-R-and-Bioconductor/Draghici/p/book/9781439809754.
Efron, Bradley. 2013. Large-scale inference : empirical Bayes methods for estimation, testing, and prediction. Cambridge University Press. http://admin.cambridge.org/academic/subjects/statistics-probability/statistical-theory-and-methods/large-scale-inference-empirical-bayes-methods-estimation-testing-and-prediction.
Gentleman, R. 2017. Annotate: Annotation for Microarrays.
Gentleman, R., V. Carey, W. Huber, and F. Hahne. 2017. Genefilter: Genefilter: Methods for Filtering Genes from High-Throughput Experiments.
Gregory Alvord, W., J. A. Roayaei, O. A. Quinones, and K. T. Schneider. 2007. “A microarray analysis for differential gene expression in the soybean genome using Bioconductor and R.” Briefings in Bioinformatics 8 (6): 415–31. https://doi.org/10.1093/bib/bbm043.
Hackstadt, Amber J, and Ann M Hess. 2009. “Filtering for increased power for microarray data analysis.” BMC Bioinformatics 10 (January): 11. https://doi.org/10.1186/1471-2105-10-11.
Huber, W., Carey, V. J., Gentleman, R., Anders, et al. 2015. “Orchestrating High-Throughput Genomic Analysis with Bioconductor.” Nature Methods 12 (2): 115–21. http://www.nature.com/nmeth/journal/v12/n2/full/nmeth.3252.html.
Irizarry, Rafael A., Bridget Hobbs, Francois Collin, Yasmin D. Beazer‐Barclay, Kristen J. Antonellis, Uwe Scherf, and Terence P. Speed. 2003. “Exploration, Normalization, and Summaries of High Density Oligonucleotide Array Probe Level Data.” Biostatistics 4 (2): 249–64. https://doi.org/10.1093/biostatistics/4.2.249.
Jeanmougin, Marine, Aurelien de Reynies, Laetitia Marisa, Caroline Paccard, Gregory Nuel, and Mickael Guedj. 2010. “Should We Abandon the t-Test in the Analysis of Gene Expression Microarray Data: A Comparison of Variance Modeling Strategies.” Edited by Kerby Shedden. PLoS ONE 5 (9): e12336. https://doi.org/10.1371/journal.pone.0012336.
Kauffmann, Audrey, Robert Gentleman, and Wolfgang Huber. 2009. “arrayQualityMetrics–a Bioconductor Package for Quality Assessment of Microarray Data.” Bioinformatics 25 (3): 415–16.
Khatri, Purvesh, and Sorin Draghici. 2005. “Ontological analysis of gene expression data: current tools, limitations, and open problems.” Bioinformatics (Oxford, England) 21 (18): 3587–95. https://doi.org/10.1093/bioinformatics/bti565.
Khatri, Purvesh, Marina Sirota, and Atul J. Butte. 2012. “Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges.” PLOS Computational Biology 8 (2): e1002375. https://doi.org/10.1371/journal.pcbi.1002375.
MacDonald, James W. 2017. Mogene21sttranscriptcluster.db: Affymetrix Mogene21 Annotation Data (Chip Mogene21sttranscriptcluster).
Sánchez-Pla, Alex. 2014. “DNA Microarrays Technology: Overview and Current Status.” In, edited by Alejandro Cifuentes Carolina Simó and Virginia García-Cañas, Volume 63:1–23. Fundamentals of Advanced Omics Technologies: From Genes to Metabolites. Elsevier. http://www.sciencedirect.com/science/article/pii/B9780444626516000015.
Slowikowski, Kamil. 2017. Ggrepel: Repulsive Text and Label Geoms for ’Ggplot2’. https://CRAN.R-project.org/package=ggrepel.
Smyth, G. K. 2005. “limma: Linear Models for Microarray Data.” In Bioinformatics and Computational Biology Solutions Using r and Bioconductor, 397–420. New York: Springer-Verlag. https://doi.org/10.1007/0-387-29362-0_23.
Smyth, Gordon K. 2004. “Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments.” Statistical Applications in Genetics and Molecular Biology 3 (1): 1–25. https://doi.org/10.2202/1544-6115.1027.
Tan, Yuande, and Yin Liu. 2011. “Comparison of methods for identifying differentially expressed genes across multiple conditions from microarray data.” Bioinformation 7 (8): 400–404. http://www.ncbi.nlm.nih.gov/pubmed/22347782 http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC3280440.
Tusher, V G, R Tibshirani, and G Chu. 2001. “Significance analysis of microarrays applied to the ionizing radiation response.” Proceedings of the National Academy of Sciences of the United States of America 98 (9): 5116–21. https://doi.org/10.1073/pnas.091062498.
Warnes, Gregory R., Ben Bolker, Lodewijk Bonebakker, Robert Gentleman, Wolfgang Huber Andy Liaw, Thomas Lumley, Martin Maechler, et al. 2016. Gplots: Various r Programming Tools for Plotting Data. https://CRAN.R-project.org/package=gplots.
Wickham, Hadley. 2009. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. http://ggplot2.org.